#pragma once

#include <vector>
#include <algorithm>
#include <cmath>
#include <memory>

namespace entropys {

inline int intpow(int x, int e)
{
    if (e == 0) return 1;
    if (x == 0) return 0;
    int r = x;
    for (int i = 1; i < e; i++)
        r *= x;
    return r;
}

inline int doublecmp(const void *a, const void *b)
{
    if (*(double*)a < *(double*)b)
        return -1;
    else
        return 1;
}

inline double sde(const double X[], int N, int epsilon, int m, int lambda)
{
    double* x = new double[N];
    memcpy(x, X, N * sizeof(double));
    qsort(x, N, sizeof(double), doublecmp);
    double* partings = new double[epsilon - 1];
    for (int i = 0; i < epsilon - 1; i++)
        partings[i] = x[N / epsilon * (i + 1) - 1];
    delete[] x;
    int* symbols = new int[N]();
    for (int i = 0; i < N; i++) {
        for (int j = epsilon - 2; j >= 0; j--) {
            if (X[i] >= partings[j]) {
                symbols[i] = j + 1;
                break;
            }
        }
    }
    delete[] partings;
    int permutations = intpow(epsilon, m);
    int* counts = new int[permutations]();
    int* counts2 = new int[permutations]();
    double* statemode = new double[permutations]();
    int** count_transition = new int* [permutations];
    double** state_transition = new double* [permutations];
    for (int i = 0; i < permutations; i++) {
        count_transition[i] = new int[epsilon]();
        state_transition[i] = new double[epsilon];
    }
    int subN = N - (m - 1) * lambda;
    for (int i = 0; i < subN; i++) {
        int index = 0;
        int j;
        for (j = 0; j < m; j++)
            index += symbols[i + j * lambda] * intpow(epsilon, m - j - 1);
        counts[index]++;
        if (i < subN - lambda) {
            counts2[index]++;
            count_transition[index][symbols[i + 1 + (j - 1) * lambda]]++;
        }
    }
    for (int i = 0; i < permutations; i++) {
        statemode[i] = 1.0 * counts[i] / subN;
        for (int j = 0; j < epsilon; j++) {
            if (count_transition[i][j] != 0)
                state_transition[i][j] = 1.0 * count_transition[i][j] / counts2[i];
            else
                state_transition[i][j] = 0;
        }
    }
    delete[] symbols;
    delete[] counts;
    delete[] counts2;
    double result = 0;
    for (int i = 0; i < permutations; i++) {
        if (statemode[i] != 0)
            result -= statemode[i] * log(statemode[i]);
        for (int j = 0; j < epsilon; j++) {
            if (state_transition[i][j] != 0)
                result -= statemode[i] * state_transition[i][j] * log(state_transition[i][j]);
        }
    }
    delete[] statemode;
    for (int i = 0; i < permutations; i++) {
        delete[] count_transition[i];
        delete[] state_transition[i];
    }
    delete[] count_transition;
    delete[] state_transition;
    return result / log(intpow(epsilon, m + 1));
}

inline double msde(const double X[], int N, int tau, int epsilon, int m, int lambda)
{
    int y_len = N - tau + 1;
    double* y = new double[y_len];
    for (int j = 0; j < y_len; j++) {
        double sum = 0;
        for (int i = j; i < j + tau; i++) {
            sum += X[i];
        }
        y[j] = sum / tau;
    }
    double result = sde(y, y_len, epsilon, m, lambda);
    delete[] y;
    return result;
}

inline double sde(const std::vector<double> &X, int epsilon, int m, int lambda)
{
    return sde(X.data(), X.size(), epsilon, m, lambda);
}

inline double msde(const std::vector<double> &X, int tau, int epsilon, int m, int lambda)
{
    return msde(X.data(), X.size(), tau, epsilon, m, lambda);
}

}